Compute area and thickness changes on a grid (experimental!)

from oggm import cfg
from oggm import tasks, utils, workflow, graphics
import xarray as xr
import matplotlib.pyplot as plt

This notebook proposes a method for redistributing glacier ice that has been simulated along the flowline after a glacier retreat simulation. Extrapolating the glacier ice onto a map involves certain assumptions and trade-offs. Depending on the purpose, different choices may be preferred. For example, higher resolution may be desired for visualization compared to using the output for a hydrological model. It is possible to add different options to the final function to allow users to select the option that best suits their needs.

This notebook demonstrates the redistribution process using a single glacier. Its purpose is to initiate further discussion before incorporating it into the main OGGM code base (currently in the sandbox).

Pick a glacier

# Initialize OGGM and set up the default run parameters
cfg.initialize(logging_level='WARNING')

# Local working directory (where OGGM will write its output)
# WORKING_DIR = utils.gettempdir('OGGM_distr4')
cfg.PATHS['working_dir'] = utils.get_temp_dir('OGGM_distributed', reset=True)

2023-06-22 09:15:20: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-06-22 09:15:20: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-06-22 09:15:20: oggm.cfg: Multiprocessing: using all available processors (N=2)

rgi_ids = ['RGI60-11.01450']  # This is Aletsch
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/elev_bands/W5E5'

gdir = workflow.init_glacier_directories(rgi_ids, prepro_base_url=base_url, from_prepro_level=3, prepro_border=80)[0]

2023-06-22 09:15:22: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2023-06-22 09:15:22: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

Experiment: currently only a random warming simulation

There is a little bit more work needed on OGGM for more complex workflows including past and future simulations (not much! coming soon).

# Do a random run with a bit of warming
tasks.run_random_climate(gdir, nyears=100, 
                         y0=2000,  # Climate of 1985-2015
                         seed=1,  # Change for another randomness 
                         temperature_bias=1,  # casual warming - change for other scenarios
                         store_fl_diagnostics=True,  # important! This will be needed for the redistribution
                         output_filesuffix='_rdn_1',  # optional - here I just want to make things explicit as to which run we are using afterwards
                        );

Redistribute: preprocessing

The required tasks can be found in the distribute_2d module of the sandbox:

from oggm.sandbox import distribute_2d

# This is to add a new topography to the file (smoothed differently)
distribute_2d.add_smoothed_glacier_topo(gdir)
# This is to get the bed map at the start of the simulation
tasks.distribute_thickness_per_altitude(gdir)
# This is to prepare the glacier directory for the interpolation (needs to be done only once)
distribute_2d.assign_points_to_band(gdir)

Let’s have a look at what we just did:

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()

# Inititial glacier thickness
f, ax = plt.subplots()
ds.distributed_thickness.plot(ax=ax);
ax.axis('equal');

# Which points belongs to which band, and then within one band which are the first to melt
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ds.band_index.plot(ax=ax1);
ds.rank_per_band.plot(ax=ax2);
ax1.axis('equal'); ax2.axis('equal'); plt.tight_layout();

Redistribute simulation

The tasks above need to be run only once. The next one however should be done for each simulation:

distribute_2d.distribute_thickness_from_simulation(gdir, input_filesuffix='_rdn_1')

Plot

Let’s have a look!

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()

ds.time

<xarray.DataArray 'time' (time: 101)>
array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
        12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,
        24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,
        36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,
        48.,  49.,  50.,  51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,
        60.,  61.,  62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,
        72.,  73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,
        84.,  85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,
        96.,  97.,  98.,  99., 100.])
Coordinates:
  * time     (time) float64 0.0 1.0 2.0 3.0 4.0 ... 96.0 97.0 98.0 99.0 100.0
Attributes:
    description:  Floating hydrological year

xarray.DataArray

'time'

• time: 101

• 0.0 1.0 2.0 3.0 4.0 5.0 6.0 ... 94.0 95.0 96.0 97.0 98.0 99.0 100.0

  array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
          12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,
          24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,
          36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,
          48.,  49.,  50.,  51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,
          60.,  61.,  62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,
          72.,  73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,
          84.,  85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,
          96.,  97.,  98.,  99., 100.])

• Coordinates: (1)
  • time
    (time)
    float64
    0.0 1.0 2.0 3.0 ... 98.0 99.0 100.0

    description :

    Floating hydrological year

    array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
            12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,
            24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,
            36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,
            48.,  49.,  50.,  51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,
            60.,  61.,  62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,
            72.,  73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,
            84.,  85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,
            96.,  97.,  98.,  99., 100.])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(Index([  0.0,   1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,
           ...
            91.0,  92.0,  93.0,  94.0,  95.0,  96.0,  97.0,  98.0,  99.0, 100.0],
          dtype='float64', name='time', length=101))

• Attributes: (1)

  description :

  Floating hydrological year

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
ds.simulation_distributed_thickness_rdn_1.sel(time=0).plot(ax=ax1, vmax=400);
ds.simulation_distributed_thickness_rdn_1.sel(time=40).plot(ax=ax2, vmax=400);
ds.simulation_distributed_thickness_rdn_1.sel(time=80).plot(ax=ax3, vmax=400);
ax1.axis('equal'); ax2.axis('equal'); plt.tight_layout();

Animation!

from matplotlib import animation
from IPython.display import HTML, display

# Get a handle on the figure and the axes
fig, ax = plt.subplots()

thk = ds['simulation_distributed_thickness_rdn_1']

# Plot the initial frame. 
cax = thk.isel(time=0).plot(ax=ax,
    add_colorbar=True,
    cmap='viridis',
    vmin=0, vmax=350,
    cbar_kwargs={
        'extend':'neither'
    }
)
ax.axis('equal')

def animate(frame):
    ax.set_title(f'Year {int(frame)}')
    cax.set_array(thk.values[frame, :].flatten())

ani_glacier = animation.FuncAnimation(fig, animate, frames=len(thk.time), interval=200);

HTML(ani_glacier.to_jshtml())

Once Loop Reflect

# Write to mp4?
# FFwriter = animation.FFMpegWriter(fps=10)
# ani2.save('animation.mp4', writer=FFwriter)